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An expansion method for perturbation of the zero 
temperature grand canonical density matrix is intro- 
duced. The method achieves quadratically conver- 
gent recursions that yield the response of the zero 
temperature density matrix upon variation of the 
Hamiltonian. The technique allows treatment of em- 
bedded quantum subsystems with a computational 
cost scaling linearly with the size of the perturbed 
region, 0(N pel t.), and as 0(1) with the total system 
size. It also allows direct computation of the den- 
sity matrix response functions to any order with lin- 
ear scaling effort. Energy expressions to 4th order 
based on only first and second order density matrix 
response are given. 

In electronic structure theory, significant effort has 
been devoted to the development of methods with the 
computational cost scaling linearly with system size 
[1,2]. The ability to perform accurate calculations with 
reduced-complexity O(N) scaling is an important break- 
through that opens a variety of new possibilities in com- 
putational materials science, chemistry and biology. One 
of the most elegant and efficient approaches to linear 
scaling is density matrix purification [3-9], where con- 
structing the density matrix by quadratically convergent 
spectral projections replaces the single-particle eigen- 
value problem arising in tight-binding and self-consistent 
Hartree-Fock and Kohn-Sham theory. For large insulat- 
ing systems this method is efficient because of the sparse 
O(N) real-space matrix representation of operators. In- 
stead of cubic scaling, the computational cost scales lin- 
early with the system size. Apart from O(N) purifica- 
tion techniques there are alternative approaches such as 
constrained functional minimization [10,11], and hybrid 
schemes [12-16]. 

In this letter, we introduce a grand canonical den- 
sity matrix perturbation theory based on recently de- 
veloped spectral projection methods for purification of 
the density matrix [8,9]. The method provides direct so- 
lution of the zero temperature density matrix response 
upon variation of the Hamiltonian through quadrati- 
cally convergent recursions. The method makes it pos- 
sible to study embedded quantum subsystems and den- 
sity matrix response functions within linear scaling effort. 
The density matrix perturbation technique avoids using 



wavefunction formalism. In spirit, it is therefore simi- 
lar to the density matrix perturbation method proposed 
by McWeeny [17] and offers a flexibility comparable to 
Green's function methods [18,19]. The present work is 
likewise related to the recent work of Bowler and Gillan 
[20], who developed a functionally constrained density 
matrix minimization scheme for embedding. However, 
our approach to computation of the density matrix re- 
sponse is quite different from existing methods of solu- 
tions for the coupled-perturbed self-consistent-field equa- 
tions. In contrast to conventional methods that pose so- 
lution implicitly through coupled equations [21-24], the 
new method provides explicit construction of the deriva- 
tive density matrix through recursion. 

The main problem in constructing a density matrix 
perturbation theory is the non-analytic relation between 
the zero temperature density matrix and the Hamilto- 
nian, given by the discontinuous step function [25], 



P = %/ - H] 



(1) 



This discontinuity makes expansion of P about H dif- 
ficult. At finite temperatures the discontinuity disap- 
pears and we may use perturbation expansions of the 
analytic Fermi-Dirac distribution [26]. However, even at 
finite temperatures a perturbation expansion based on 
the Fermi-Dirac distribution may have slow convergence. 

In linear scaling purification schemes [3-9] , the density 
matrix is constructed by recursion; 



X 
X n 

p 



= L(H), 

-1 = Fn(X n ) 

= linu,^ 



n = 0,1,2,. 



(2) 



Here, L{H) is a linear normalization function mapping 
all eigenvalues of H in reverse order to the interval of con- 
vergence [0, 1] and F n (X n ) is a set of functions projecting 
the eigenvalues of X n toward 1 (for occupied states) or 
(for unoccupied states) . In one of the simplest and most 
efficient techniques, which requires only knowledge of the 
number of occupied states N e and no a priori knowledge 
of n [8], we have 



P n (X n ) 



Y 2 

IX, 



n 



Tr(X n ) > N e 
Tr{X n ) < N e . 



(3) 
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Purification expansion schemes are quadratically conver- 
gent, numerically stable, and can even solve problems 
with degenerate eigenstates and fractional occupancy [9] . 
Thanks to an exponential decay of the density matrix el- 
ements as a function of |r — r' for insulating materials, 
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the operators have a sparse matrix representation and 
the number of non-zero matrix elements above a numer- 
ical threshold scales linearly with the system size. In 
these cases the matrix-matrix multiplications, which are 
the most time consuming steps, have an TV-scaling cost. 

Equivalent to the purification schemes are the sign- 
matrix expansions [27-29]. The general scheme is the 
same as in Eq. (2), but the expansion is performed around 
a step from — 1 to 1 at x = 0. 

Our grand canonical density matrix perturbation the- 
ory is based on the purification in Eq. (2). A perturbed 
Hamiltonian H = + gives the expansion 



X n =X^+A n , n = 0,l,2, 



(4) 



where X n ^ is the unperturbed expansion and A„ are the 
differences due to the perturbation It is then easy 

to show that 



— F n (Xn^ + A„) — F n (Xn^) 

P = p(°)+lim rwoc A„. 



(5) 



This is the key result of the present article and defines 
our grand canonical density matrix perturbation theory. 
Combining Eq. (5) with the expansion in Eq. (3) gives 
the recursive expansion [25] 



{X n 0) ,A n } + A 2 n if Tr(Xk u ')>N t 
2A„ - {X n °\ A„} - A 2 n otherwise. 
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A n 
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(6) 



Other expansions based on, for example, McWeeny, trace 
conserving or trace resetting purification [3,5,9] can also 
be included in this quite general approach. However, 
Eq. (6) is particularly efficient since it only requires two 
matrix multiplications per iteration. Because the pertur- 
bation expansions inherit properties from their generator 
sequence, they are likewise quadratically convergent with 
iteration, numerically stable, and exact to within accu- 
racy of the drop tolerance [9] . 

If the perturbed xjf ' has eigenvalues outside the in- 
terval of convergence [0, 1] the expansion could fail. To 
avoid this problem the normalization function L(H) in 
Eq. (2) can be chosen to contract the eigenvalues of Xq 
to [S, 1 — 5} , where S > is sufficiently large. 

A major advantage with the expansion in Eq. (6) is 
that for band-gap materials that are locally perturbed, 
the A„ are likewise localized as a result of nearsighted- 
ness [30,11]. The matrix products in Eq. (6) can therefore 
be calculated using only the local regions of X n that re- 
spond to the perturbation. Given that perturbation does 
not change the overall decay of the density matrix, the 
computational cost of the expansion scales linearly with 
the size of the perturbed region 0(N pelt .) and as 0(1) 
with the total system size. 

Density matrix purification docs not necessarily re- 
quire prior knowledge of the chemical potential, but once 



the initial expansion of the unperturbed system is carried 
out, the chemical potential is set. The perturbation ex- 
pansions of Eq. (5) are therefore grand canonical [31]. 
With this in mind, Eq. (6) may be readily applied to em- 
bedding schemes that do involve long range charge flow. 

The computation of many spectroscopic properties 
such as the Raman spectra, chemical shielding and polar- 
ization requires the calculation of density matrix deriva- 
tives with respect to perturbation. Grand canonical den- 
sity matrix perturbation theory can be used to compute 
these response functions. Assume a perturbation of the 
Hamiltonian H^°\ 



h = h<v + \hW, 



(7) 



in the limit A — > 0. The corresponding density matrix is 
P = P (°) +\ P W +\2 P W + ... , (g) 

where the response functions pW (density matrix deriva- 
tives) correspond to order /x in A. Expanding the pertur- 
bation as in Eq. (6), individual response terms may be 
collected order by order at each iteration; 



A„ = AA« + \ 2 A^ + 



(9) 



Keeping terms through order m in A at each iteration, 
with A^ = X n , the following recursive sequence is ob- 
tained for fjb = m, m — 1, . . . , 1 : 



" +1 1 2A^-£f =0 A«A<r 4) 



if Tr(X n ) > N e 



otherwise. 



(10) 



These equations provide an explicit, quadratically con- 
vergent solution of the response functions, where 



= lim A'"'. 



(11) 



With the same technique it is possible to treat perturba- 



r(2) 



tions where H = ff(°> + \ a m i} + \ b H^ } + \ a \ b H { a 
to produce a mixed density matrix expansion P = p(°) + 



A Q Pi x) + \ b P?> + X a \ b P% 

Equation (10) provide direct explicit construction of 
the response equations based on well developed linear 
scaling technologies [8,9]. This is quite different from 
conventional approaches [21-24] that pose solution im- 
plicitly through coupled matrix equations, achieving at 
best linear scaling with iterative solvers. 

Higher order expansions of an observable can be cal- 
culated efficiently from low order density matrix terms. 
Similar to Wigner's 2n+l rule for wavefunctions [32] we 
have the energy response E = P (0) + AP (1) + \ 2 E (2) + 
X 3 E^ + \ A E^\ where 

= Tr(P^H^), E<® = 0.5Tr(P«P«) 
E^ =Tr([pW,P<V]pWHV>), (12) 
B< 4 ) - 0.5Tr([(2I - p(°))p( 2 )p(°)p(i) 

_p(0)p(l)p(2)(j + p(0) )]#(!)). 



(1) 



,(2) 
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The corresponding n+1 rule for /j, > is 



5 n +i = S 1 dS(S 1 + <*>„), 



(14) 



E W = M - 1 Tr(P ( ^- 1) i? (1) ). (13) 

To demonstrate the grand canonical density matrix 
perturbation theory, we present two examples based on 
single-site perturbations of a model Hamiltonian and a 
beta-carotene molecule. 

The model Hamiltonian has random diagonal elements 
exhibiting exponential decay of the overlap elements as a 
function of site separation on a randomly distorted lat- 
tice. This model represents a Hamiltonian of an insu- 
lator that might occur, for example, with a Gaussian 
basis set in density functional theory or in various tight- 
binding schemes. A local perturbation is imposed on 
the model Hamiltonian by moving the position of one 
of the lattice sites. Using the perturbation expansion of 
Eq. (6), a series of perturbations A„ is generated. In 
each step a numerical threshold t = 10~ 6 is applied 
as described in [9]. The lower inset in Fig. 1 shows 
the number of elements above the threshold in A„ as 
a function of iteration. The local perturbation is effi- 
ciently represented with only ~ 50 elements out of 10 4 . 
Figure 1 also illustrates the quadratic convergence of 
the error. At convergence after M iterations the new 
perturbed density matrix is given by P = p(°) + Am- 
The error \\P — P 0X act||2 = 6.4 x 1CP 5 and the error 
\E - Poxactl = 1.3 x 1(T 6 [34]. The error of the per- 
turbed density matrix P is stable at convergence and 
close to the numerical error for the unperturbed density 
matrix due to thresholding | |p(°)-P ( c °j act 1 1 2 = l.OxlCT 4 , 
and|p(°)-££L t | = 2.4xl(r 6 . 

The electronic structure of the second example, the 
beta-carotene molecule, was calculated with the Mon- 
doSCF suite of linear scaling algorithms [33] at the 
RHF/STO-2G level of theory. Figure 2 shows the matrix 
sparsity factor of the density matrix for beta-carotene. 
The difference between two fully self-consistent Fockians 
was chosen as a perturbation, (one with and one with- 
out a small displacement of a single carbon atom). In 
this way, more long-ranged effects due to self-consistency 
are included. Even if beta-carotene is too small to have 
a very sparse representation of the density matrix, the 
perturbation sequence A„ is found to be highly sparse. 
The error with threshold r = 10 -5 in the single-particle 
Hartree-Fock energy \E — P cxa ct| = 2.8 x 1CP 5 a.u., which 
is of the same order of error as for the unperturbed 
molecule. Standard first order perturbation theory yields 
an error two orders of magnitude larger [35]. 

The present formulation has been developed in an or- 
thogonal representation. With a TV-scaling congruence 
transformation [12], it is straightforward to employ this 
representation when using a non-orthogonal basis. When 
using a non-orthogonal basis set, change in the inverse 
overlap matrix S^ 1 due to a local perturbation dS is 
given by the recursive Dyson equation, 



where S = S - dS, 5 = 0, and S' 1 = So' 1 + 
lim^oo 5 n . The equation contains only terms with local 
sparse updates and the computational cost scales linearly, 
0(N pert .), with the size of the perturbed region. Similar 
perturbation schemes for the sparse inverse Cholesky or 
square root factorizations can also be applied [36] . 

Density matrix perturbation theory can be applied in 
many contexts. For example, a straightforward calcu- 
lation of the energy difference due to a small pertur- 
bation of a very large system may not be possible be- 
cause of the numerical problem in resolving a tiny en- 
ergy difference between two large energies. With density 
matrix perturbation theory, we work directly with the 
density matrix difference A„ and the problem can be 
avoided, for example, the single particle energy change 
AP = linin^oo Tr(HA n ). In analogy to incremen- 
tal Fock builds in self-consistent field calculations [39], 
the technique can be used in incremental density ma- 
trix builds. Connecting and disconnecting individual 
weakly interacting [40] quantum subsystems can be per- 
formed by treating off-diagonal elements of the Hamil- 
tonian as a perturbation. This should be highly use- 
ful in nanoscience for connecting quantum dots, surfaces, 
clusters and nanowires, where the different parts can be 
calculated separately, provided a connection through a 
common chemical potential is given, for example via a 
surface substrate. In quantum molecular dynamics, such 
as quantum mechanical-molecular mechanical QM/MM 
schemes, or Monte-Carlo simulations, where only a local 
part of the system is perturbed and updated, the new ap- 
proach is of interest. Several techniques used within the 
Green's function context also should apply for the den- 
sity matrix. The proposed perturbation approach may be 
used for response functions, impurities, effective medium 
and local scattering techniques [18,19,37,38]. The theory 
of grand canonical density matrix perturbation is thus 
a rich field with applications in many areas of materials 
science, chemistry and biology. 

In summary, we have introduced a grand canonical 
perturbation theory for the zero temperature density 
matrix, extending quadratically convergent purification 
techniques to expansions of the density matrix upon vari- 
ation of the Hamiltonian. The perturbation method al- 
lows the local adjustment of embedded quantum sub- 
systems with a computational cost that scales as 0(1) 
for the total system size and as 0(N pcrt .) for the re- 
gion that respond to the perturbation, as demonstrated 
in Figs. 1 and 2. A new quadratically convergent N- 
scaling recursive approach to computing density matrix 
response functions has been outlined, and energy expres- 
sions to 4th order in terms of only first and second order 
density matrix response were given. The density matrix 
perturbation technique is surprisingly simple and offers 
an efficient alternative to several Green's function meth- 
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ods and conventional schemes for solution of the coupled 
perturbed self-consistent-field equations. 

Discussions with E. Chisolm, S. Corish, S. Tretiak, C. 
J. Tymczak, V. Weber, and J. Wills are gratefully ac- 
knowledged. 



Corresponding author: amn@lanl.gov 

S. Goedecker, Rev. Mod. Phys. 71, 1085 (1999). 

S. Y. Wu and C. S. Jayanthi, Physics Reports-Review 

Section of Physics Letters, 358, 1 (2002). 

R. McWeeny, Rev. Mod. Phys. 32, 335 (I960). 

W. L. Clinton, A. J. Galli, and L. J. Masa, Phys. Rev. 

177, 7 (1969). 

A. H. R. Palser and D. E. Manolopoulos, Phys. Rev. B 
58, 12704 (1998). 

A. Holas, Chem. Phys. Lett, 340, 552 (2001). 

A. M. N. Niklasson, C. J. Tymczak, and H. Roder, Phys. 

Rev. B 66, 155120 (2002). 

A. M. N. Niklasson, Phys. Rev. B 66, 155115 (2002). 
A. M. N. Niklasson, C. J. Tymczak, and M. Challacombe, 
J. Chem. Phys. 118, 8611 (2003). 

X. -P. Li, W. Nunes, and D. Vanderbilt, Phys. Rev. B 
47, 10891 (1993). 

W. Kohn, Phys. Rev. Lett. 76, 3168 (1996). 

M. Challacombe, J. Chem. Phys. 110, 2332 (1999). 

D. R. Bowler and M. J. Gillan, Comput. Phys. Commun. 

120, 95 (1999). 

A. D. Daniels and G. E. Scuseria, J. Chem. Phys. 110, 
1321 (1999). 

T. Helgaker, H. Larsen, J. Olsen, P Jorgensen, Chem. 
Phys. Lett. 327, 397 (2000). 

M. Head-Gordon, Y. H. Shao, C. Saravanan, C. A. 
White, Mol. Phys. 101, 37 (2003); Y. H. Shao, C. Sar- 
avanan, M. Head-Gordon, C. A. White, J. Chem. Phys. 
118, 6144 (2003). 

R. McWeeny, Phys. Rev. 126, 1028 (1962). 

R. Haydock, Solid State Phys. Advances in Research and 

Applications, 35, 215 (1980). 

J. E. Inglesfield, J. Phys. C, 14, 3795 (1981). 

D. R. Bowler and M. J. Gillan, Chem. Phys. Lett. 355, 

306 (2002). 

M. Frisch, M. Head-Gordon, and J. Pople, Chem. Phys. 
141, 306 (1990). 

S. P. Kama and M. Dupuis, J. Comput. Chem. 12, 487 
(1991). 

C. Ochsenfeld and M. Head-Gordon, Chem. Phys. Lett. 
270, 399 (1997). 

H. Larsen, P. Jorgensen, and T. Helgaker, J. Chem. Phys. 
113, 8908 (2000). 

The notation: 9 is the Heaviside step function, /i 
is the chemical potential, I — S(r — r'), H — 
H(r,r'), P = P(r,r'), Tr(A) = JdrA(r,r), AB = 
J dr"A(r,r")B(r",r'), the density n(r) = P(r,r), the 
single particle energy E = Tr(HP), N e is the number 
of occupied states below /i, and N is the total num- 



ber of states corresponding to the system size. The anti- 
commutator {A, B} — AB + BA and the commutator 
[A, B] = AB — BA. 

[26] R. P. Feynam, "Statistical Mechanics", ISBN 0-201- 
36076-4, Addison Wesley Longman, Inc. (1998). 

[27] C. S. Kenney, and A. J. Laub, SIAM J. Matrix Anal. 2, 
273 (1991). 

[28] G. Beylkin, N. Coult, and M. Mohlenkamp, J. Comput. 

Phys. 152, 32 (1999). 
[29] K. Nemeth, and G. S. Scuseria, J. Chem. Phys. 113, 6035 

(2000) . 

[30] W. Kohn, Phys. Rev. 115, 809 (1959). 

[31] A canonical perturbation theory can be constructed by 
allowing the spectral projections F n (X n + A„) to differ 
from F n (X„) in Eq. (5), however, in this case the locality 
of the expansion for a local perturbation is lost. 

[32] T. Helgaker, P. Jorgensen, and J. Olsen, Molecular 
electronic- structure theory, ISBN 0471 96755 6, (John 
Wiley & Sons, West Sussex, England, 2000). 

[33] M. Challacombe, E. Schwengler, C. J. Tymczak, C. K. 
Gan, K. Nemeth, V. Weber and A. M. N. Niklasson, 
MondoSCF A program suite for massively parallel, lin- 
ear scaling SCF theory and ab initio molecular dynamics 

(2001) , URL http://www.tl2.lanl.gov/home/mchalla/. 
[34] As a comparison, first-order perturbation theory gives an 

error \Tr((H (0) + AV)P (0) ) - £Wact = 8.7 x 10~ 3 . 
[35] |Tr((F (0) + AV)P (0) ) - Bexact] = 4.5 x 10" 3 a.u., where 

F^ is the unperturbed Fockian. 
[36] A. M. N. Niklasson, (unpublished). 

[37] I. Turek, V. Drachl, J. Kudrnovsky, M. Sob, and P. Wein- 
berger, Electronic structure of disordered alloys, surfaces 
and interfaces, Klawer Academic Publishers (1997). 

[38] I. A. Abrikosov, A. M. N. Niklasson, S. I. Simak, B. Jo- 
hansson, A. V. Ruban, and H. L. Skriver, Phys. Rev. 
Lett. 76, 4203 (1996). 

[39] E. Schwengler, M. Challacombe, and M. Head-Gordon, 
J. Chem. Phys. 106, 9708 (1997). 

[40] In the case of strong interaction and self-consistent cal- 
culations, long ranged coupling may occur from charge 
fluctuations and changed dipole field. In this case the 
perturbation is no longer local. 



FIG. 1. The Error = log 10 (||Al 0) + A n - P eX act|| 2 ) as a 
function of iterations n (N = 100, N e = 50). The lower 
inset shows the number of non-zero matrix elements of A„ 
above threshold r = 10~ 6 . The upper inset shows the density 
matrix perturbation. 



FIG. 2. The matrix sparsity factors (number of non-zero 
elements over the total number of elements) for beta-carotene 
(RHF/STO-2G). 
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